Isomap — Nonlinear Dimensionality Reduction (Geodesic Distances)

Isomap (Isometric Mapping) is a manifold learning method that tries to preserve walking distances along a surface, not straight-line “as the crow flies” distances.

What you’ll learn

• why Euclidean distance can be misleading on curved manifolds (Swiss roll)

• how Isomap works: kNN graph → geodesic distances → classical MDS

• how to implement Isomap from scratch (NumPy + shortest paths)

• how to diagnose common failure modes (noise, disconnected graphs)

import heapq

import numpy as np
import plotly.express as px
import plotly.graph_objects as go
import os
import plotly.io as pio

from plotly.subplots import make_subplots
from scipy.linalg import orthogonal_procrustes
from sklearn.manifold import Isomap as SklearnIsomap
from sklearn.metrics import pairwise_distances
from sklearn.neighbors import NearestNeighbors


pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")
np.set_printoptions(precision=4, suppress=True)

rng = np.random.default_rng(7)

Prerequisites

• Euclidean distance + basic linear algebra

• Graphs and shortest paths (Dijkstra / Floyd–Warshall)

• Eigendecomposition (classical MDS)

1) Intuition: “Measuring distance by walking paths, not flying”

If you’re in a city with rivers and bridges, two places can be:

• close “as the crow flies” (straight line)

• but far if you must walk along streets (path distance)

Isomap is built on the same idea:

• Euclidean distance ≈ “flying”

• geodesic distance ≈ “walking along the manifold”

Swiss-roll metaphor

A Swiss roll is a 2D sheet rolled into 3D.

• Points on different layers can be close in 3D Euclidean distance.

• But along the sheet, they can be far.

Isomap tries to “unroll” the sheet by preserving geodesic distances.

def make_swiss_roll(n_samples: int = 350, noise: float = 0.0, random_state: int = 7):
    rng_local = np.random.default_rng(random_state)

    # Roughly matches sklearn's swiss-roll parameterization
    t = 1.5 * np.pi * (1.0 + 2.0 * rng_local.random(n_samples))  # [1.5π, 4.5π]
    h = 21.0 * rng_local.random(n_samples)

    x = t * np.cos(t)
    z = t * np.sin(t)
    y = h
    X = np.column_stack([x, y, z])

    if noise > 0:
        X = X + rng_local.normal(scale=noise, size=X.shape)

    return X, t, h


X, t, h = make_swiss_roll(n_samples=350, noise=0.0, random_state=7)
X.shape

(350, 3)

fig = px.scatter_3d(
    x=X[:, 0],
    y=X[:, 1],
    z=X[:, 2],
    color=t,
    title="Swiss roll in 3D (color = intrinsic coordinate t)",
    labels={"x": "x", "y": "y (height)", "z": "z", "color": "t"},
)
fig.update_traces(marker=dict(size=4, opacity=0.85))
fig.update_layout(scene=dict(aspectmode="data"))
fig.show()

fig = px.scatter(
    x=t,
    y=h,
    color=t,
    title="Ground truth 2D coordinates (t, height)",
    labels={"x": "t (along the sheet)", "y": "height", "color": "t"},
)
fig.update_traces(marker=dict(size=6, opacity=0.85))
fig.show()

2) Graph-based explanation

Isomap is a graph + geometry method:

1. Build a k-nearest neighbors (kNN) graph in the original space.

2. Use shortest paths in that graph to approximate geodesic distances.

3. Run classical MDS on those geodesic distances to get a low-dimensional embedding.

2.1 kNN graph

• each point is a node

• connect each node to its k nearest neighbors (by Euclidean distance)

• edge weight = that local Euclidean distance

Locally, Euclidean distance is a good approximation to distance along the manifold.

def knn_graph(X: np.ndarray, k: int) -> list[list[tuple[int, float]]]:
    """Undirected kNN graph with edge weights = Euclidean distances."""
    nn = NearestNeighbors(n_neighbors=k + 1, algorithm="auto")
    nn.fit(X)
    dists, idx = nn.kneighbors(X)

    idx = idx[:, 1:]
    dists = dists[:, 1:]

    n = X.shape[0]
    nbr_dicts: list[dict[int, float]] = [dict() for _ in range(n)]
    for i in range(n):
        for j, d in zip(idx[i], dists[i]):
            j = int(j)
            w = float(d)
            nbr_dicts[i][j] = min(nbr_dicts[i].get(j, np.inf), w)
            nbr_dicts[j][i] = min(nbr_dicts[j].get(i, np.inf), w)

    return [sorted(nbrs.items(), key=lambda t: t[1]) for nbrs in nbr_dicts]


def connected_components(neighbors: list[list[tuple[int, float]]]):
    n = len(neighbors)
    comp = -np.ones(n, dtype=int)
    cid = 0

    for start in range(n):
        if comp[start] != -1:
            continue
        stack = [start]
        comp[start] = cid
        while stack:
            u = stack.pop()
            for v, _ in neighbors[u]:
                if comp[v] == -1:
                    comp[v] = cid
                    stack.append(v)
        cid += 1

    return comp, cid


def ensure_connected_knn_graph(X: np.ndarray, k_start: int = 4, k_max: int = 30):
    for k in range(k_start, k_max + 1):
        neighbors = knn_graph(X, k=k)
        _, n_comp = connected_components(neighbors)
        if n_comp == 1:
            return neighbors, k
    raise ValueError(f"Graph stayed disconnected for k in [{k_start}, {k_max}]")


k_target = 10
neighbors, k = ensure_connected_knn_graph(X, k_start=k_target, k_max=30)
k

10

# Visualize the kNN graph on a smaller subset (edges get cluttered fast)
n_vis = 160
idx_vis = rng.choice(X.shape[0], size=n_vis, replace=False)

Xv = X[idx_vis]
tv = t[idx_vis]

neighbors_v, k_v = ensure_connected_knn_graph(Xv, k_start=10, k_max=25)

xs, ys, zs = [], [], []
for a, nbrs in enumerate(neighbors_v):
    for b, _w in nbrs:
        if a < b:
            xs.extend([Xv[a, 0], Xv[b, 0], None])
            ys.extend([Xv[a, 1], Xv[b, 1], None])
            zs.extend([Xv[a, 2], Xv[b, 2], None])

fig = go.Figure()
fig.add_trace(
    go.Scatter3d(
        x=xs,
        y=ys,
        z=zs,
        mode="lines",
        line=dict(color="rgba(0,0,0,0.18)", width=1),
        name=f"kNN edges (k={k_v})",
        hoverinfo="skip",
    )
)
fig.add_trace(
    go.Scatter3d(
        x=Xv[:, 0],
        y=Xv[:, 1],
        z=Xv[:, 2],
        mode="markers",
        marker=dict(size=4, color=tv, colorscale="Viridis", opacity=0.9, colorbar=dict(title="t")),
        name="points",
    )
)
fig.update_layout(
    title=f"kNN graph on Swiss roll (n={n_vis}, k={k_v})",
    scene=dict(aspectmode="data"),
)
fig.show()

2.2 Geodesic distances via shortest paths

Once we have a neighborhood graph, we approximate geodesic distances by computing shortest path distances in the graph:

• Use Dijkstra from every source node (common when edges are non-negative).

• Or use Floyd–Warshall (simple, but \(O(n^3)\); good for small graphs).

These shortest-path distances are our estimate of “walking distance” along the manifold.

def dijkstra_distances(neighbors: list[list[tuple[int, float]]], source: int):
    n = len(neighbors)
    dist = np.full(n, np.inf)
    dist[source] = 0.0

    heap: list[tuple[float, int]] = [(0.0, source)]
    while heap:
        d_u, u = heapq.heappop(heap)
        if d_u != dist[u]:
            continue
        for v, w_uv in neighbors[u]:
            nd = d_u + w_uv
            if nd < dist[v]:
                dist[v] = nd
                heapq.heappush(heap, (nd, v))

    return dist


def dijkstra_with_prev(neighbors: list[list[tuple[int, float]]], source: int):
    n = len(neighbors)
    dist = np.full(n, np.inf)
    prev = -np.ones(n, dtype=int)

    dist[source] = 0.0
    heap: list[tuple[float, int]] = [(0.0, source)]
    while heap:
        d_u, u = heapq.heappop(heap)
        if d_u != dist[u]:
            continue
        for v, w_uv in neighbors[u]:
            nd = d_u + w_uv
            if nd < dist[v]:
                dist[v] = nd
                prev[v] = u
                heapq.heappush(heap, (nd, v))

    return dist, prev


def all_pairs_shortest_paths(neighbors: list[list[tuple[int, float]]]):
    n = len(neighbors)
    D = np.empty((n, n), dtype=float)
    for s in range(n):
        D[s] = dijkstra_distances(neighbors, source=s)
    return D


def adjacency_matrix(neighbors: list[list[tuple[int, float]]]):
    n = len(neighbors)
    W = np.full((n, n), np.inf)
    np.fill_diagonal(W, 0.0)
    for i, nbrs in enumerate(neighbors):
        for j, w in nbrs:
            if w < W[i, j]:
                W[i, j] = w
    return W


def floyd_warshall(W: np.ndarray):
    D = W.copy()
    n = D.shape[0]
    for k in range(n):
        D = np.minimum(D, D[:, [k]] + D[[k], :])
    return D

D_geo = all_pairs_shortest_paths(neighbors)
D_euc = pairwise_distances(X)

np.isfinite(D_geo).all(), D_geo.shape

(True, (350, 350))

# A quick sanity check that Floyd–Warshall matches Dijkstra (on a smaller graph)
idx_small = rng.choice(X.shape[0], size=120, replace=False)
X_small = X[idx_small]
neighbors_small, k_small = ensure_connected_knn_graph(X_small, k_start=10, k_max=25)

D_dij = all_pairs_shortest_paths(neighbors_small)
W_small = adjacency_matrix(neighbors_small)
D_fw = floyd_warshall(W_small)

np.max(np.abs(D_dij - D_fw))

7.105427357601002e-15

3) Plotly visualizations

3.1 Euclidean vs geodesic distances (matrix view)

If we sort points by the intrinsic coordinate \(t\), a Swiss roll has a clear “rolled up” structure.

• Euclidean distances can connect across layers.

• Geodesic distances tend to respect the manifold.

order = np.argsort(t)
m = 120
idx_h = order[:m]

E = D_euc[np.ix_(idx_h, idx_h)]
G = D_geo[np.ix_(idx_h, idx_h)]

fig = make_subplots(
    rows=1,
    cols=2,
    subplot_titles=("Euclidean distance (flying)", "Geodesic distance (walking)"),
)
fig.add_trace(go.Heatmap(z=E, colorscale="Viridis", showscale=False), row=1, col=1)
fig.add_trace(go.Heatmap(z=G, colorscale="Viridis", showscale=True), row=1, col=2)
fig.update_layout(title=f"Euclidean vs geodesic distances (subset m={m}, sorted by t)")
fig.show()

3.2 Two points: “close if you fly”, “far if you walk”

We’ll automatically find two points that have:

• small Euclidean distance, but

• large difference in intrinsic coordinate \(t\) (likely different layers)

Then we compare:

• the straight-line chord (Euclidean)

• the shortest path along the kNN graph (geodesic estimate)

def find_cross_layer_pair(D_euc: np.ndarray, t: np.ndarray, min_delta_t: float = 7.0):
    D = D_euc.copy()
    np.fill_diagonal(D, np.inf)

    for thr in [min_delta_t, 6.0, 5.0]:
        mask = np.abs(t[:, None] - t[None, :]) > thr
        cand = np.where(mask, D, np.inf)
        i, j = np.unravel_index(np.argmin(cand), cand.shape)
        if np.isfinite(cand[i, j]):
            return int(i), int(j), float(thr)

    raise ValueError("No suitable cross-layer pair found")


i, j, thr = find_cross_layer_pair(D_euc, t, min_delta_t=7.0)

dist_i, prev_i = dijkstra_with_prev(neighbors, source=i)

path = [j]
cur = j
while cur != i and cur != -1:
    cur = int(prev_i[cur])
    path.append(cur)
path = path[::-1]

print(f"Chosen pair: i={i}, j={j} (|Δt| > {thr})")
print(f"Euclidean distance: {D_euc[i, j]:.3f}")
print(f"Geodesic distance (graph shortest path): {D_geo[i, j]:.3f}")
print(f"Number of edges in path: {len(path) - 1}")

Chosen pair: i=176, j=196 (|Δt| > 7.0)
Euclidean distance: 8.904
Geodesic distance (graph shortest path): 25.947
Number of edges in path: 5

path_xyz = X[path]

fig = go.Figure()
fig.add_trace(
    go.Scatter3d(
        x=X[:, 0],
        y=X[:, 1],
        z=X[:, 2],
        mode="markers",
        marker=dict(size=3, color=t, colorscale="Viridis", opacity=0.55),
        name="points",
    )
)
fig.add_trace(
    go.Scatter3d(
        x=[X[i, 0], X[j, 0]],
        y=[X[i, 1], X[j, 1]],
        z=[X[i, 2], X[j, 2]],
        mode="lines",
        line=dict(color="red", width=7),
        name="Euclidean chord",
    )
)
fig.add_trace(
    go.Scatter3d(
        x=path_xyz[:, 0],
        y=path_xyz[:, 1],
        z=path_xyz[:, 2],
        mode="lines+markers",
        line=dict(color="black", width=4),
        marker=dict(size=4, color="black"),
        name="Graph shortest path",
    )
)
fig.add_trace(
    go.Scatter3d(
        x=[X[i, 0]],
        y=[X[i, 1]],
        z=[X[i, 2]],
        mode="markers",
        marker=dict(size=9, color="red"),
        name="start",
    )
)
fig.add_trace(
    go.Scatter3d(
        x=[X[j, 0]],
        y=[X[j, 1]],
        z=[X[j, 2]],
        mode="markers",
        marker=dict(size=9, color="red"),
        name="end",
    )
)
fig.update_layout(
    title=f"Close in Euclidean, far in geodesic (k={k})",
    scene=dict(aspectmode="data"),
)
fig.show()

4) Classical MDS on geodesic distances (the “Isomap” embedding step)

Classical MDS takes a distance matrix \(D\) and finds points in a low-dimensional Euclidean space whose pairwise distances match \(D\) as closely as possible.

Given distances \(D_{ij}\), classical MDS:

1. squares distances: \(D^2\)

2. double-centers them to build an inner-product (Gram) matrix \(B\)

3. takes the top eigenvectors/eigenvalues of \(B\)

Isomap is simply:

classical MDS, but using geodesic distances instead of Euclidean ones.

Implementation note: to match the supervised-learning notebooks, we’ll wrap the full pipeline in a small ScratchIsomap class (fit → geodesic distances → MDS embedding).

def classical_mds(D: np.ndarray, n_components: int = 2):
    if not np.isfinite(D).all():
        raise ValueError("Distance matrix has non-finite entries (graph disconnected?)")

    n = D.shape[0]
    D2 = D**2

    J = np.eye(n) - np.ones((n, n)) / n
    B = -0.5 * J @ D2 @ J

    eigvals, eigvecs = np.linalg.eigh(B)
    idx = np.argsort(eigvals)[::-1]
    eigvals = eigvals[idx]
    eigvecs = eigvecs[:, idx]

    eigvals = np.maximum(eigvals, 0.0)
    Y = eigvecs[:, :n_components] * np.sqrt(eigvals[:n_components])
    return Y, eigvals


class ScratchIsomap:
    def __init__(
        self,
        n_neighbors: int = 10,
        n_components: int = 2,
        shortest_path: str = "dijkstra",
        ensure_connected: bool = True,
        k_max: int = 30,
    ):
        self.n_neighbors = int(n_neighbors)
        self.n_components = int(n_components)
        self.shortest_path = str(shortest_path)
        self.ensure_connected = bool(ensure_connected)
        self.k_max = int(k_max)

    def fit(self, X: np.ndarray):
        X = np.asarray(X, dtype=float)
        if X.ndim != 2:
            raise ValueError('X must be a 2D array')
        n = X.shape[0]
        if n < 2:
            raise ValueError('X must have at least 2 samples')
        if self.n_components < 1:
            raise ValueError('n_components must be >= 1')
        if self.n_neighbors < 1:
            raise ValueError('n_neighbors must be >= 1')
        if self.n_neighbors >= n:
            raise ValueError('n_neighbors must be < n_samples')

        if self.ensure_connected:
            neighbors, k_used = ensure_connected_knn_graph(X, k_start=self.n_neighbors, k_max=self.k_max)
        else:
            neighbors = knn_graph(X, k=self.n_neighbors)
            k_used = self.n_neighbors

        if self.shortest_path == 'dijkstra':
            D_geo = all_pairs_shortest_paths(neighbors)
        elif self.shortest_path in {'floyd-warshall', 'floyd_warshall'}:
            D_geo = floyd_warshall(adjacency_matrix(neighbors))
        else:
            raise ValueError("shortest_path must be 'dijkstra' or 'floyd_warshall'")

        Y, eigvals = classical_mds(D_geo, n_components=self.n_components)

        self.X_ = X
        self.neighbors_ = neighbors
        self.n_neighbors_ = k_used
        self.geodesic_distances_ = D_geo
        self.embedding_ = Y
        self.eigvals_ = eigvals
        return self

    def fit_transform(self, X: np.ndarray) -> np.ndarray:
        return self.fit(X).embedding_


scratch_iso = ScratchIsomap(n_neighbors=k, n_components=2, shortest_path='dijkstra', ensure_connected=False).fit(X)

# Use the scratch model outputs downstream
neighbors = scratch_iso.neighbors_
k = scratch_iso.n_neighbors_
D_geo = scratch_iso.geodesic_distances_
Y = scratch_iso.embedding_
eigvals = scratch_iso.eigvals_

eigvals[:10]

array([42278.283 , 38156.4554, 23052.3745, 16412.0213,  4571.1708,
        4122.0284,  2861.4723,  2175.4576,  1825.6325,  1413.0827])

fig = px.line(
    x=np.arange(1, 11),
    y=eigvals[:10],
    markers=True,
    title="Top eigenvalues (classical MDS on geodesic distances)",
    labels={"x": "component", "y": "eigenvalue"},
)
fig.show()

fig = make_subplots(rows=1, cols=2, subplot_titles=("Ground truth (t, h)", "Isomap embedding (from scratch)"))

fig.add_trace(
    go.Scatter(
        x=t,
        y=h,
        mode="markers",
        marker=dict(size=4, color=t, colorscale="Viridis", opacity=0.8),
        showlegend=False,
    ),
    row=1,
    col=1,
)

fig.add_trace(
    go.Scatter(
        x=Y[:, 0],
        y=Y[:, 1],
        mode="markers",
        marker=dict(
            size=4,
            color=t,
            colorscale="Viridis",
            opacity=0.8,
            colorbar=dict(title="t"),
        ),
        showlegend=False,
    ),
    row=1,
    col=2,
)

fig.update_xaxes(title_text="t", row=1, col=1)
fig.update_yaxes(title_text="h (height)", row=1, col=1)
fig.update_xaxes(title_text="component 1", row=1, col=2)
fig.update_yaxes(title_text="component 2", row=1, col=2)

fig.update_layout(title="Unrolling the Swiss roll: Isomap vs ground truth")
fig.show()

5) Algorithm steps (summary)

1. Neighborhood graph

   • build a kNN graph (nodes = points, weighted edges = local Euclidean distances)

2. Shortest paths

   • compute all-pairs shortest paths to approximate geodesic distances

   • common choices: Dijkstra (run from each node) or Floyd–Warshall (small graphs)

3. Embedding

   • run classical MDS on the geodesic distance matrix to get low-dimensional coordinates

5.1 Swiss-roll unfolding animation (Plotly)

To make the “unrolling” feel concrete, we’ll animate a transition from:

• the original 3D Swiss roll, to

• a 2D Isomap embedding placed onto a flat plane

This is just a visual interpolation, not part of the algorithm.

def zscore(A: np.ndarray):
    return (A - A.mean(axis=0)) / A.std(axis=0)


# Align the Isomap embedding to the "true" unrolled coordinates for a nicer animation.
T_true = np.column_stack([t, h])
Y0 = zscore(Y)
T0 = zscore(T_true)

R, _ = orthogonal_procrustes(Y0, T0)
Y_aligned = Y0 @ R

Xc = X - X.mean(axis=0)
scale = float(np.mean(Xc.std(axis=0)))
plane = np.column_stack([Y_aligned[:, 0], Y_aligned[:, 1], np.zeros(X.shape[0])]) * scale

mins = np.minimum(Xc.min(axis=0), plane.min(axis=0))
maxs = np.maximum(Xc.max(axis=0), plane.max(axis=0))
pad = 0.06 * (maxs - mins)
ranges = [(mins[i] - pad[i], maxs[i] + pad[i]) for i in range(3)]

alphas = np.linspace(0.0, 1.0, 31)

def frame_points(alpha: float):
    P = (1.0 - alpha) * Xc + alpha * plane
    return P


P0 = frame_points(alphas[0])
trace0 = go.Scatter3d(
    x=P0[:, 0],
    y=P0[:, 1],
    z=P0[:, 2],
    mode="markers",
    marker=dict(size=4, color=t, colorscale="Viridis", opacity=0.9, cmin=t.min(), cmax=t.max(), colorbar=dict(title="t")),
)

frames = []
for a in alphas:
    P = frame_points(float(a))
    frames.append(
        go.Frame(
            data=[
                go.Scatter3d(
                    x=P[:, 0],
                    y=P[:, 1],
                    z=P[:, 2],
                    mode="markers",
                    marker=dict(size=4, color=t, colorscale="Viridis", opacity=0.9, cmin=t.min(), cmax=t.max()),
                )
            ],
            name=f"{a:.2f}",
        )
    )

steps = [
    {
        "args": [[fr.name], {"frame": {"duration": 0, "redraw": True}, "mode": "immediate"}],
        "label": fr.name,
        "method": "animate",
    }
    for fr in frames
]

fig = go.Figure(data=[trace0], frames=frames)
fig.update_layout(
    title="Swiss-roll unfolding (interpolating to a 2D Isomap embedding plane)",
    scene=dict(
        xaxis=dict(range=ranges[0], title="x"),
        yaxis=dict(range=ranges[1], title="y"),
        zaxis=dict(range=ranges[2], title="z"),
        aspectmode="data",
    ),
    updatemenus=[
        {
            "type": "buttons",
            "showactive": False,
            "x": 0.05,
            "y": 1.08,
            "buttons": [
                {
                    "label": "Play",
                    "method": "animate",
                    "args": [
                        None,
                        {
                            "frame": {"duration": 80, "redraw": True},
                            "fromcurrent": True,
                            "transition": {"duration": 0},
                        },
                    ],
                },
                {
                    "label": "Pause",
                    "method": "animate",
                    "args": [[None], {"frame": {"duration": 0, "redraw": False}, "mode": "immediate"}],
                },
            ],
        }
    ],
    sliders=[
        {
            "active": 0,
            "currentvalue": {"prefix": "unfold α="},
            "pad": {"t": 45},
            "steps": steps,
        }
    ],
)
fig.show()

6) Failure modes

6.1 Sensitivity to noise

Isomap relies on the kNN graph approximating the true manifold locally.

• With too much noise, local neighborhoods get corrupted.

• That makes geodesic distances (shortest paths) unreliable.

6.2 Disconnected graphs

If k is too small, the kNN graph can become disconnected.

• Some geodesic distances become infinite.

• Classical MDS can’t embed an infinite distance matrix.

# Noise example
Xn, tn, hn = make_swiss_roll(n_samples=350, noise=1.0, random_state=7)
neighbors_n, k_n = ensure_connected_knn_graph(Xn, k_start=10, k_max=30)
D_geo_n = all_pairs_shortest_paths(neighbors_n)
Yn, _ = classical_mds(D_geo_n, n_components=2)

fig = make_subplots(rows=1, cols=2, subplot_titles=("Isomap (clean)", "Isomap (noise=1.0)"))

fig.add_trace(
    go.Scatter(
        x=Y[:, 0],
        y=Y[:, 1],
        mode="markers",
        marker=dict(size=4, color=t, colorscale="Viridis", opacity=0.8),
        showlegend=False,
    ),
    row=1,
    col=1,
)

fig.add_trace(
    go.Scatter(
        x=Yn[:, 0],
        y=Yn[:, 1],
        mode="markers",
        marker=dict(size=4, color=tn, colorscale="Viridis", opacity=0.8, colorbar=dict(title="t")),
        showlegend=False,
    ),
    row=1,
    col=2,
)

fig.update_xaxes(title_text="component 1", row=1, col=1)
fig.update_yaxes(title_text="component 2", row=1, col=1)
fig.update_xaxes(title_text="component 1", row=1, col=2)
fig.update_yaxes(title_text="component 2", row=1, col=2)

fig.update_layout(title="Noise can distort neighborhoods and geodesic distances")
fig.show()

# Disconnected graph example (k too small)
neighbors_bad = knn_graph(X, k=3)
comp_bad, n_comp_bad = connected_components(neighbors_bad)
sizes = np.bincount(comp_bad)

print(f"Connected components: {n_comp_bad}")
print("Component sizes:", sizes)

fig = px.scatter_3d(
    x=X[:, 0],
    y=X[:, 1],
    z=X[:, 2],
    color=comp_bad.astype(str),
    title="kNN graph can disconnect when k is too small (here: k=3)",
    labels={"x": "x", "y": "y", "z": "z", "color": "component"},
)
fig.update_traces(marker=dict(size=4, opacity=0.85))
fig.update_layout(scene=dict(aspectmode="data"))
fig.show()

D_geo_bad = all_pairs_shortest_paths(neighbors_bad)
finite_fraction = float(np.isfinite(D_geo_bad).mean())
print(f"Finite geodesic distance fraction: {finite_fraction:.3f}")

if not np.isfinite(D_geo_bad).all():
    print("Disconnected graph → some geodesic distances are infinite. Increase k (or use a radius graph).")

Connected components: 2
Component sizes: [344   6]

Finite geodesic distance fraction: 0.966
Disconnected graph → some geodesic distances are infinite. Increase k (or use a radius graph).

7) Practical usage: scikit-learn’s Isomap

scikit-learn’s Isomap does the same conceptual steps (kNN graph → shortest paths → embedding), but with a well-tested implementation.

iso = SklearnIsomap(n_neighbors=k, n_components=2)
Y_sklearn = iso.fit_transform(X)

fig = make_subplots(rows=1, cols=2, subplot_titles=("From scratch", "sklearn"))
fig.add_trace(
    go.Scatter(
        x=Y[:, 0],
        y=Y[:, 1],
        mode="markers",
        marker=dict(size=4, color=t, colorscale="Viridis", opacity=0.8),
        showlegend=False,
    ),
    row=1,
    col=1,
)
fig.add_trace(
    go.Scatter(
        x=Y_sklearn[:, 0],
        y=Y_sklearn[:, 1],
        mode="markers",
        marker=dict(size=4, color=t, colorscale="Viridis", opacity=0.8, colorbar=dict(title="t")),
        showlegend=False,
    ),
    row=1,
    col=2,
)

fig.update_xaxes(title_text="component 1", row=1, col=1)
fig.update_yaxes(title_text="component 2", row=1, col=1)
fig.update_xaxes(title_text="component 1", row=1, col=2)
fig.update_yaxes(title_text="component 2", row=1, col=2)
fig.update_layout(title=f"Isomap embeddings (k={k})")
fig.show()

Exercises

1. Try different values of k. When does the embedding start to “short-circuit” (connect across layers)?

2. Increase noise gradually. How much noise can Isomap tolerate before it fails to unroll the Swiss roll?

3. Replace the kNN graph with an ε-radius graph. How does that change connectivity and the embedding?

References

• Tenenbaum, de Silva, Langford (2000): A Global Geometric Framework for Nonlinear Dimensionality Reduction

• scikit-learn docs: sklearn.manifold.Isomap